Codecov Report

❌ Patch coverage is 48.14815% with 28 lines in your changes missing coverage. Please review.
✅ Project coverage is 97.04%. Comparing base (bc1048b) to head (1cf5115).

@@            Coverage Diff             @@
##             main    #2910      +/-   ##
==========================================
- Coverage   97.09%   97.04%   -0.05%     
==========================================
  Files         630      631       +1     
  Lines       48855    48898      +43     
==========================================
+ Hits        47435    47450      +15     
- Misses       1420     1448      +28     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

We need jlchan/StartUpDG.jl#219 to be merged and released before.

Finally, we probably still want to figure out why elixir_advection_restart_amr.jl has a slightly different result. I can try to find that out.

I bisected this discrepancy to this PR: SciML/OrdinaryDiffEq.jl#3441, which first landed in OrdinaryDiffEqCore.jl v3.33.0. Using the commit directly before that PR:

julia> include("examples/tree_2d_dgsem/elixir_advection_restart_amr.jl");
[...]
julia> analysis_callback(sol)
(l2 = [8.018497923389368e-5], linf = [0.0007307237754662355])

(run_v3_32) pkg> st
Status `~/.julia/dev/Trixi/run_v3_32/Project.toml`
  [7d9f7c33] Accessors v0.1.44
  [bbf590c4] OrdinaryDiffEqCore v3.32.0 `https://github.com/SciML/OrdinaryDiffEq.jl.git:lib/OrdinaryDiffEqCore#97f41902`
⌃ [b0944070] OrdinaryDiffEqLowStorageRK v1.15.0
  [a7f1ee26] Trixi v0.16.8-DEV `..`
Info Packages marked with ⌃ have new versions available and may be upgradable.

With the commit from the PR:

julia> include("examples/tree_2d_dgsem/elixir_advection_restart_amr.jl");
[...]
julia> analysis_callback(sol)
(l2 = [8.018498574373939e-5], linf = [0.0007307237754662355])

(run_v6) pkg> st
Status `~/.julia/dev/Trixi/run_v6/Project.toml`
  [7d9f7c33] Accessors v0.1.44
  [bbf590c4] OrdinaryDiffEqCore v3.32.0 `https://github.com/SciML/OrdinaryDiffEq.jl.git:lib/OrdinaryDiffEqCore#1864296`
⌃ [b0944070] OrdinaryDiffEqLowStorageRK v1.15.0
  [a7f1ee26] Trixi v0.16.8-DEV `..`
Info Packages marked with ⌃ have new versions available and may be upgradable.

These L2-error values are exactly the two values in https://github.com/trixi-framework/Trixi.jl/pull/2910/changes#diff-0d913e7961860bd3a9c595d68dbd7a7b09ea17758afc8c814b9e45a454a04796R77. So this has nothing to do with restarting, controllers, or the v7 updates of the SciML packages, but is simply a coincidence. It didn't surface before since we currently restrict OrdinaryDiffEqCore.jl to v3.31. Looking at the changes of that PR, I think we can simply use the new value. Do you agree, @ranocha?

Thanks for debugging this! I agree that we can modify the value in this case.

So the new getters in SciML/OrdinaryDiffEq.jl#3570 should resolve these or the cache setup step?

With SciML/OrdinaryDiffEq.jl#3570 being merged and released, the correct default values are picked up again and I reverted the changes related to that (and bumped the compat bound for OrdinaryDiffEqCore.jl accordingly). So I think, we only need to wait for #3003 to be merged until this PR is ready. Let me know if you need help with that PR @vchuravy. I would like to finish this, such that we can start integrating the new versions into downstream packages like TrixiShallowWater.jl and TrixiAtmo.jl.

Failing CI is due to changes in OrdinaryDiffEqSDIRK.jl. With OrdinaryDiffEqSDIRK.jl v2.1.0 I can reproduce the reference values and with v2.2.0 I get similar results as CI, which don't pass the tolerance. It's very likely this comes from SciML/OrdinaryDiffEq.jl#3196. Two concrete examples are:

• Running trixi_include(joinpath(examples_dir(), "tree_1d_dgsem", "elixir_advection_diffusion.jl"), tspan=(0.0, 0.4)); (using KenCarp4(autodiff = AutoFiniteDiff(), linsolve = LUFactorization())) gives with OrdinaryDiffEqSDIRK.jl v2.2.0

julia> analysis_callback(sol)
(l2 = [8.389498175038025e-6], linf = [2.8474220044816256e-5])

and with OrdinaryDiffEqSDIRK.jl v2.1.0

julia> analysis_callback(sol)
(l2 = [8.40483031802723e-6], linf = [2.8990878868540015e-5])

The difference is bigger than the tolerances of the ode integrator (absolute and relative tolerances of 1e-10):

julia> 8.389498175038025e-6 - 8.40483031802723e-6
-1.533214298920447e-8

julia> 2.8474220044816256e-5 - 2.8990878868540015e-5
-5.16658823723759e-7

• Running

trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem",
                                        "elixir_navierstokes_viscous_shock_newton_krylov.jl"),
                     tspan=(0.0, 0.1),
                     atol_lin_solve=1e-11,
                     rtol_lin_solve=1e-11,
                     atol_ode_solve=1e-10,
                     rtol_ode_solve=1e-10);

(using linsolve = KrylovJL_GMRES(atol = atol_lin_solve, rtol = rtol_lin_solve) ; ode_alg = Kvaerno4(autodiff = AutoFiniteDiff(), linsolve = linsolve)) gives with OrdinaryDiffEqSDIRK.jl v2.2.0

julia> analysis_callback(sol)
(l2 = [0.005262765896012976, 0.00544494709949278, 2.3365757290728324e-17, 0.0055117947406083685], linf = [0.010910620789879921, 0.011624105849751709, 3.1032813029796676e-16, 0.011845859748702448])

and with OrdinaryDiffEqSDIRK.jl v2.1.0

julia> analysis_callback(sol)
(l2 = [3.428501006913763e-5, 2.5967418005918904e-5, 2.7662844381999672e-17, 2.855861765164612e-5], linf = [0.00018762342908829055, 0.00014059002077504434, 3.7252310624872113e-16, 0.00014510700486725092])

which have the following differences:

julia> analysis_callback(sol).l2 .- l2
4-element StaticArraysCore.SizedVector{4, Float64, Vector{Float64}} with indices SOneTo(4):
  0.005228480885943838
  0.005418979681486861
 -4.297087091271349e-18
  0.005483236122956722

julia> analysis_callback(sol).linf .- linf
4-element StaticArraysCore.SizedVector{4, Float64, Vector{Float64}} with indices SOneTo(4):
  0.01072299736079163
  0.011483515828976665
 -6.219497595075436e-17
  0.011700752743835197

It's also worth noting that OrdinaryDiffEqSDIRK.jl v2.2.0 takes much smaller time steps (~7.9728e-05) than v2.1.0 (~1.0726e-03).
Can you help, @ChrisRackauckas?

@singhharsh1708 can you look into this?

Is the difference only apparent with Krylov linear solvers? If so, that would narrow down what is possibly the issue by a lot.

No, the first example uses LUFactorization() for linsolve.

@JoshuaLampert Tried to reproduce on a minimal heat-equation MWE with KenCarp4 + AutoFiniteDiff + LUFactorization (with and without a mass matrix) at abstol=reltol=1e-10. Pre-#3196 master, current master, and the open @generated PR (#3642) all give matching step counts and l2 errors (differences only at FP precision). Could you share a smaller (non-Trixi) repro, or re-check against v2.3.0 / #3642?

Thanks, @singhharsh1708! I will retry with v2.3.0 and SciML/OrdinaryDiffEq.jl#3643. Creating a non-Trixi.jl MWE could be hard and could take some time.

With v2.3.0 the situation already improved and we only have one more failing example (the second I posted above). I also tested it locally with SciML/OrdinaryDiffEq.jl#3643 and got the same as with v2.3.0. Could you maybe look for some issue specific to Kvaerno4 or KrylovJL_GMRES or the combination of both @singhharsh1708?

@JoshuaLampert Still working on it — I’ll keep you posted.

I locally verified that the remaining failing parabolic tests are fixed by SciML/OrdinaryDiffEq.jl#3653.

I think moving forward parts of the IO logic (for restarts) should be moved into a sciml package, as

1. any PDE package will benefit from having restartable integrators
2. maintenance of the IO can be simplified this way, because failures are discovered earlier and so are changes in potential IO invariants